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SIMIARY 

The behaviour of the steady state spurious error modes of the MacComiack scheme 
and the ipwind scheme of Warming and Beam is obtained from a linearized difference 
equation for the steady state error. It is shown that the spurious errors can 
exist either as an eigensolution of the hanogeneous part of this difference 
equation or because of excitation from large discretization errors near oblique 
shocks. It is found that the upwind scheme does not permit spurious oscillations 
on the upstream side of shocks. Exanples are given for the inviscid Burgers' 
equation and for one- and two-dimensional gasdynamic flows. 
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1. INTRODUaiO^ 


Che of the major difficulties in numerical conpitaticHi of gasdynamic flows is 
caused by the occurrence of shock waves. The most straightforward way of dealing 
with shocks is based on the concept of weak solutions of conservation laws in- 
troduced by Lax tH • By a conservative formulation of the difference equations 
the shocks are captured automatically and no special treatment is required. The 
price to be paid for this simplicity is generally large errors in the vicinity of 
the shodcs, most commonly in the form of pre- and post-shock oscillations. Effi- 
cient use of these so called shock-capturing methods therefore requires som under- 
standing of the nature of these errors. 

The existing work on this question is mostly based on the time-dependent equations. 
One standard procedure is to investigate the truncation error of simple linear 
model equations [2]. The results obtained are limited to the amplification and 
phase errors of the Fourier components of the solution and do not give any 
direct information on the effects of shocks. The case of an unsteady travelling 
shock has been treated by Kreiss and Lundqvist [3] for the simple linear wave 
equation 3u/ 3t - 3u/ 3x * 0 and by Lerat and Peyret [4] for the nonlinear, 
inviscid Burgers’ equation. Both these works show that oscillating errors are 
generally to be expected near the shock. 

Time-dependent numerical methods are also a useful tool for many steady flows of 
aerodynamic interest. In this paper the errors in the final steady state of such 
computations will be considered. Since none of the above methods of analysis are 
particularly well suited for the steady case, we shall describe a procedure by 
vhich it is possible to study the steady state errors more directly [7] , [14] • 
Consider the system of conservation laws in one space dimension 


3F ^ 3H 

5T ^ 


0 


( 1 ) 


where for inviscid, compressible flow F= {p,m, e } is the field vector of mass, 
x-momentum, and total energy densities. We solve this system by the stable time- 
-dependent numerical scheme 

n^l 

F. = Dj (F^'') j = 1, 2, ... , N (2) 

where Fj is the computed field vector at grid point j at time step n , 
and D j is the difference operator replacing the partial differential problem , 
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i.e. eq. (1) with appropriate boundary conditions. The procedure (2) has con- 
verged to a steady state if ^ ^ ^ specified accuracy in all the N 

grid points. Our purpose is to study the error n. this result, particularly when 
shocks appear in the field. The error is defined by: 


e*? » - f". 


(3) 


vdiere F^. = F (n At, j Ax) is the exact solution of the differential problem with 
• • • 0 

initial conditions Fqj. Substitution of Fj from eq. (3) into the difference 
(2) gives an equation for the error expressed in the exact solution Fqj . 

Such an equation is of course just as conplicated as the original difference pro- 
blem and furthermore requires a knowledge of the unknown exact solution Fq* The 
discussion is therefore limited to small errors, and the nonlinear effects will 
be neglected. Substitution of (3) into (2) gives after linearization: 



L. (e^) 

3 ^ J 



(4) 


viiere Lj is the linearized operator and 


T? = Dj (i^j) - 


■Oj 


(5) 


is the truncation error. From (4) the usual requirements for convergence of the 
linearized problem is obtained. If the error shall become smaller as the mesh is 
refined, the operator must be consistent so that 0. This is particularly 

important if shocks appear in the field. We must also require the operator Lj to 
be stable [2]. When a steady state is reached, eqs. (4) and (5) reduce to: 


(I - L. ) e! = T? 
J J J 


1 ■ “j - "oj 


( 6 ) 

(7) 


If we assume that the exact solution F^j is known, then Tj and Lj are detemined, 
and eq. (6) is a linear algebraic system for the error at every grid point. The 
scheme is consistent with p-th order of accuracy if Tj = 0 (Ax^), where Ax is 
the grid spacing. From eq. (7) it follows that the scheme must be consistent 
everywhere if the error shall vanish as Ax-^ 0. In particular eq. (7) implies 
that the difference scheme must be a consistent approximation of the junp relations 
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vhen shocks or other discontinuities appear in the field. Finally we must also 
require that the coefficient matrix obtained from the operator I-Lj remains non- 
singular as Ax 0. If this is not the case, there can be an error in the final 
steady state (of indetermined anplitude in the linear approximation) determined 
by eq. (6) vhich is a singular hranogeneous system in this limit. It is interest- 
ing to note that in this case an eigaivalue of Lj approaches 1 , so that the time- 
-dependent scheme ceases to converge vhen 0 . Below we shall give sin?)le 

exanples of such behaviour for a model equation and for steady one-aimensional 
flow with a shock. 

Eqs. (6) and (7) will be used below to study the behaviour of small errors for 
some sin 5 )le steady cases with known piecewise constant exact solutions . 

The analysis for more general Fq would be prohibitively difficult even if such 
a solution was known. For a slowly varying field these single results could be 
used to judge the errors locally. The analysis will be applied to the two-step 
scheme suggested by MacCormack [5], and also with the upwind corrector of Warming 
and Beam [ 6 ] , although some of the results are applicable to a wider class of 
schemes. 

We shall first discuss the homogeneous error problem which arises from a single 
nonlinear scalar equation. The properties of the operator I - Lj of eq. (6) will 
be discussed and an exan^le given which does not converge in the limit Ax 0 , 
although the computation is stable. An analysis of the initial value pi-oblem for 
this scalar equation has been given by Marten et al [14], for the Lax-Wendroff 
scheme. Their analysis shows that the steady state is weakly unstable in the L 2 
norm. It can be shown, however, that the initial-boundary value problem is stable 
when the number of grid points N is bounded which is also indicated by the com- 
putational results in [14]. 

Similar results will be derived for the one-dimensional gasdynamic case with a 
stationary shock. For the two-dimensional case we shall finally discuss the 
magnitude of the discretization error for steady, oblique shocks. 


2. A SCALAR MODEL EQUATION 

As an illustration we shall study the nonlinear, inviscid Burgers' equation 
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subject to the following initial and boundary conditions: 


u(x, 0) 


1 

- 0 

-1 


X « 0 

0 < X < 1 

x«l 


u(0,t) - 1 
u(l, t) --1 


The exact solution of this problem is 


Up(X,t) 


1 - S(x- ^ t) - S(x-1+ jt); t< 1 
. 1 -2S(x- i) ; til 


\^re S(x) is the unit step function. This solution represents two discontinuities 
(shocks) travelling away from the boundaries with velocities | and . When they 
meet at x= 5, they form a stationary discontinuity with the junp 2. 

For the one-dimensional case given by eq. (1) with F=u and H=ju^ the MacCormack 
(MC) scheme is: 


F^^^ = F^ - ah(Hj - (8a) 

fJ^'' = ^ (F^ + ^ - ah (f^ - )) (8b) 

where h= At /Ax, = H(F*?*^), and a = l gives backward predictor - forward 

corrector differencing, and the converse for a=-1. IVhen u is positive , the 
predictor with 0=! is upwind, and the second-order upwind scheme of Warming 
and Beam (WB) is obtained if (8b) is replaced by the corrector: 


j^i+l 

j 


= - (F^ 
2 


F^+l 


h (H^ - 2 H^ 

J J- 






- 


(8c) 
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In a steady state we get from eq. (8) for the MC (with a = 1 ) 

and WB schemes respectively: 


Hj -Hj., - 0 

(9 a) 


(9b) 


idiere the siperscripts have been dropped for sinplicity. 


TOe only possible exact steady solution in the one-dimensional case is 
Hq*H(Fq) = const. Since this solution exactly satisfies the difference eqs. (8), 
the steady state truncation error t| of eq. (7) is exactly zero. The eq. (6) 
for small errors on the steady state is therefore homogeneous. For the present 
scalar problem we substitute Uj = Ugj + into the difference eq. (9 a) : 

t ^ - .1 ) “Oj+1 ^ 

- (1 + huQj) %-l ^;-1 ' ° 

This equation, together with the boundary conditions u!|= = 1, “ "1 » 

i.e. =£j^ = 0, leads to a system of linear algebraic equations corresponding to 
eq. (6). If we substitute the exact steady solution Uq into these, we could 
investigate under which circumstances the coefficient matrix becomes singular. 

This method, however, is not very practical for more complicated situations and 
a different approach is therefore taken. Let us assume that the shock is located 
between mesh points m and m+1, and hence separates two regions of constant Ug. 
With Ugj =Ug = const. * 0 eq. (10) gives: 


(1 


•v) + 2 V ^ - (1 +v) = 0 


( 11 ) 


where v= hug is the Courant number. Eq. (11) is a linear difference equation 
with constant coefficients and has solutions of the form Cj = const. • A^. Eq. (11) 
gives two solutions = 1 and A£ general solution is : 


e. = k^ 


+ k,(A,)J 


( 12 ) 
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The mode is necessary for consistency and is termed the correct mode. The 

other mode is an additional solution of the difference scheme and is therefore 
called a spurious mode. Since the variation of the error is exponential in space, 
we shall refer to this second mode as the exponential mode. For stability of the 
scheme it is required that |v | £l. X 2 is hence always negative, and the exponen- 
tial mode is therefore oscillating. For j £ m Uq is positive, X 2 — “ ^ 
exponential mode is increasing towards the shock. For j > m+1 is negative 
and we therefore define the Courant number by v^-hUg to make it positive. In 
this region the eigenvalue X 2 is which means that the errors decay with the 
same rate away frran the shock on both sides. If we use the notation X 2 l^he 

general solution in this region is: 


e 


j 


k.( 



(13) 


To study the errors through the shock discontinuity it is more suitable to con- 
sider the error in the flux H= ju^. In the linear approximation 6j =Ugj Ej 
which makes eq. (10) an equation for , and eqs. (12) and (13) remain valid 
with £j replaced by 6j. From eq. (10) it can be seen that (12) is valid 
for j < m- 1 and (13) for j ^m+1. With j=m and ’Jom"”^Om+l rela- 

tion) eq. (10) gives a relation for the error across the shock: 


*^m+l ” "^m-l 


= 0 


(14) 


With the boundary conditions 6^ =6jj = 0 we also obtain from (12) and (13): 


"^2 ^2 

( 15 ) 

1 N 

-'<4 ( x :) 

( 16 ) 


From (12) and (13) we obtain with (15) and (16) 


1 


m-1 


*^m ~ 1 m-2 ^2 ‘^m-1 

i-(tt) 




N-m 


Xo6„ 


t-(-tt) 


^ N-m-1 ''2“m+l 


(17) 
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Eqs. (14) and (17) are two equaticms for and It can clearly be seen 

that these equations beccane identical, and hence the operator I “ L j singular, when 
the factors ( 17 X 2 )^ in (17) vanish. This happens when either v -*■ 1 or N,m 
The last possibility corresponds to Ax 0. From eq. (17) it is also seen that 
when such errors arc possible, they are dominated by the exponential modes, and 
the constants and kj vanish in the above limits. The decay of the error away 
from the shock made possible by the exponential modes are therefore the iirportant 
factor. This is of inportance in more conplicated problems v/here it would be 
difficult to study the coefficient matrix from I - . 

Fig. 1 shows the result of a numerical cai^utatior. of the present problem with 
N*20 and v* .4 after 100 time-steps. The oscillating steady state error can 
clearly be seen, although the decay of the exponential modes is relatively slow 
for this low Courant nuni>er. Because of the singular nature of the error in the 
steady state, the an^litude is mainly determined by the errors introduced in the 
confutation of the discontinuous unsteady transient. To illustrate this effect we 
have also made a confutation with initial conditions given by the exact solution 
of the differential problem at the time when the discontinuities have travelled two 
mesh intervals from the boundaries. This leads to a more well behaved transient in 
the subsequent numerical cranputation. Fig. 2 shows the final steady state, and 
the error is now about 1.5 t. The ratio of the errors on both sides of the shock 
‘^m+1 ^ ^m-1 result which is in good agreement with the value 

predicted from eq. (14) . 

In [ 7 ] we have shown that the error equations can be made nonsingular sinf ly by 
switching from o » 1 for Uq > 0 to o = - 1 for Uq < 0 in the MC scheme. Confuta- 
tion with this switching gives the exact steady solution to machine accuracy as 
expected. Since the equations for constant Uq (12) and (13) are the same as 
before, the error equations are nonsingular because of the special relation across 
the shock which does not permit the oscillating exponential modes. This type of 
switching is not very effective for more complicated systems of conservation laws, 
because we then have several coexisting exponential modes as will be seen below 
for the gasdynamic case. 

A much more efficient way to make the operator I-Lj nonsingular is to use a scheme 
with a different behaviour in tlie exponential modes. We therefore consider the 
scheme with the WB ifwind corrector (8c). Proceeding in the same way as for the MC 
scheme on the steady state equation (9b) it is easily shown that the eigenvalue of 
the exponential mode for Uy > 0 is given by: 
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X 


2 


V-1 
V- S 


This scheme is stable for v 1 2 [6] . The exponential mode is oscillating for v> 1 
and monotonic for v < 1 . The interesting fact is, hwever, that the error always 
decays with increasing j towards the shock. If the MC scheme is used as before 
for Uq 1 0 , the exponential mode is everyvhere decaying with increasing j . Since 
the error is zero at j ■ 1 , it must remain so at every point towards the shock. 

To keep the scheme fully conservative the switching operator of Cb) must be used. 
When the MC scheme is used for u < 0 , a steady oscillation is still possible on 
this side of the shock 19]. If the upwind scheme is used on both sides, the exact 
steady state is obtained as expected. 


3. ONE-DINENSIONAL INVISCID. / IMPRESSIBLE FLOW 


_2 

FOr this case we have F - { p, m, e } and H * (m , ~ + p, (e +p) - } wlicrc the 

P Pi 

pressure for a perfect gas is given by the equation of state p » (y-1) (e - j ^ ) , 


y being the ratio of specific heats. The only possible exact steady solution is 
Hq- const., « 0 and the error equations are homogeneous as for the previous 
case. Hq» const, is the Rankine - Hugoniot relation, giving a possible sipersonic 
state Fq^ and a subsonic state Fq 2 such that H(Fq^ * “ ^^0 ’ 


A small error Cj on the steady state Fq gives an error in vAich is 6^ ^j» 

provided that tte Jacobian Aj of H with respect to F is given by Aj • A(Fgj ) , 
i.e. we neglect its variation with Cj . Aj must also be nonsingular which rules 
out the special cases M*0 , 11, where is the Mach nuni)er, u is the velocity 

a 

and a the speed of sound. Substitution of this into cq. (9a) and using (8a) for 
the error in the predictor, we get the following equation for the steady state 
errors of the MC scheme: 


CI-hAj,,)Aj„ •h(A.^,.ApA.e. - (LhAj) A.., e.., - 


(18) 


We shall investigate the simple case ccMisisting of a constant upstream sigxjrsonic 
flow Fq^ (j £ m), a normal shock between j *m and j ■ m*l , and downstream sub- 
sonic flow ^02 (j 1 m+1). Since Aj is constant except for the juny across the 
shock, we first study the behaviour of the error both sides of the shock by 
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solving eq. (18) for Aj - const. - A. The solutions are of the form Cj “Eq , 
irfiich leads to the eigenvaliw problem: 

f (X*-l) I - (X*-2X +l)hA) Eq - 0 

There are obviously 3 eigenvalues X ■ 1 which represent the constant correct 
mode. The problem for the remaining 3 eigenvalues for the exponential inodes is: 


so that 

1 ♦ 1 . 

H Tpr “i 

where are the eigenvalues of A given by “a (.M+1) , 0)2 and oj^^aCM-l). 
In the following it will be asstmed that M > 0. The von Neumann CMidition for 
stability is then for the MC scheme: 

V - ha(M+1)< 1 

The eigenvalues X^ can therefore be written 



For V < 1 the eigenvalues are always negative which corresponds to oscillatory 
increasing or decaying modes. The X^ mode is always decaying in the upstream 
direction and docs so more rapidly as v 1. Fig. 3 shows X 2 and Xj as functions 

of the Mach number for different v. The X^ mode also decays in the igistream 
directicwi, whereas the aj mode decays upstream for M > 1 and downstream for M< 1. 
The X^ mode is very si^ificant since it allows errors produced at the shock to 
decay in both directions away from it. The decay rates grow towards the X'^ mode as 
M -*■ 00 . As M -► 0 the decay rate of the X^ mode vanishes since X, - 1, and the 
downstream decay of the X^ mode becomes increasingly more rapid, Xj > U . 
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The eigenvectors Cq correspwding to the eigenvalues of the exponential modes are 
given by the column vector matrix: 


1 


a (M+ 1) 


a*(y M* + a^ ^ 


1 

aM 

.2 


M* 


1 \ 

afM -1J 

a^Cy NP - M j 


v^ich are the eigenvectors of the Jacobian A. They arc always linearly 
independent for y 1 . 

Since we have taken A to be constant ;ind linearized, these results arc valid for 
any Lax- Wendroff type scheme using 3 points in the spatial coordinate. It can 
also be expected that the results should be locally valid i( A is slowly vaiying, 
such as for cxanple in a slowly diverging one -dimensional duct flow. 

In the ccxistant flow on both sides of tlic shock the cxpoiv'ntial part of the error 
can be written: 


Ej • 11 A-’ K 




where A is a diagonal matrix with , as elements ;uk 1 K is ::n ;u.ii)litude 

1 , _ , .I 

vector. A more suitable form of cq. (IDJ is: 


Cj., ' i: s 


[-III 


where O I: A li 


.-1 


As before it is better to consider tlie error in II for the cqiuitions valid through 
the shock, l-rom cq. (20) we get for the upstrciim and dowi'.strcam sides ivsjx'ct i vely 


• ‘-1 «m-l 


i- - 


Vi ■ *^2 


C2) 


A rclatiai for the error across the shock is obtained from cq. (18) with j * m : 
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(I-hA^) ^hfA2^A,)6^ -lI-hA,)6^.1 -0 

which can be identically written as 

Vl * C23) 

It is easily shown that the two sides of eq. (23) are identically zero by eq. (21) 
and eq. (22) respectively. The three eqs. (21) - (23) for the errors 6^_i> 
and are hence not linearly independent, and a non-zero solution can exist 

for the error through the shock. On the subsonic side of the shock there is only 
one mode idiich decays in the downstream directirai, given by the eigenvalue 
The other two inodes cannot exist because they would produce an increasingly large 
error at the downstream boundary vdiere the boundary conditions must be satisfied. 
The error therefore has the form given by the eigenvector of the Xj mode . 
Because of the exponential decay with j , the rem ining error at the downstream 
boundary rapidly approaches zero as the mesh is reJined, thus satisfying the exact 
downstream boundary conditions. From eq. (22) the error 6^^ in the last point 
upstream of the shock must have the same form as , only larger by the factor 

( 1/ Xj )2 . In the siperscmic part all the exponential modes decay upstream, and 
the error therefore vanishes at the upstream bondary as the mesh is refined, thus 
satisfying the given inlet conditions. Since the eigenvectors of the modes are 
linearly independent, any error can decay in this manner. 

To illustrate this error behaviour a sinple numerical experiment has been made. 

Hie mesh consisted of 50 grid points, and the initial conditions were taken as the 

exact steady solution with the shock at m = 25 , except for a 10 % increase in the 

pressure at the first downstream point j = 26 . At j = 1 the inlet conditions were 

kept fixed during the conputation with p^=l kg/m , m^=600 kg/m s, e^ = 4.3»10 J/wT, 

which gives an upstream Mach number = 1.60 . The exact values were also kept at 

the downstream boundary. Tlie conputation was run for 2000 time-steps with upstream 

Courant number v, = .9. The relative change in density per tune-step was then less 
-5 

than 10 . The resultiiig steady state relative error in the flux vector near the 

shock is shown in Fig. 4 . Table I shows this conputed relative error together with 
the errors calculated from the linear theory above. To determine the anplitude of 
the theoretical modes the error in the mass flux was chosen to have the same value 
at the first downstream point j =m+1. As can be seen the agreement is very good. 

Although the maximum relative error in the flux vector is only about .61, 

corresponding errors in density aiid pressure is 2 and 4 % respectively. 
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By similar calculations it is easily shown that the e^qxHiential modes of the WB 
scheme (8 c) are given by: 



This scheme is stable for M > 1 and v < 2 . For v < 1 the error is monotonic 
since >0. The mode is oscillating for v > 1 , and the other tvro for 
sufficiently high ^bch nunber. The behaviour of X '2 and X^ is shown in Fig. 5 . 

As for previous scalar case all the eigenvalues have the interesting property 
{Xj I <1 vdiich means that the error always decays in the streamwise direction. 
Since the error upstream in the supersonic part is cancelled by the specified inlet 
conditions, the exponential modes are absent down to the shock. If the switching 
operator of [ 6 ] is used with the MC schane downstream, the X^ mode can still 
remain on the downstream side [ 9 ] . An interesting proposal to remedy this situa- 
tion is the splitting technique of Warming and Beam [15 ] , which treats the waves 
going upstream in subsonic parts by downstream differencing similar to the WB 
scheme . 


4. TWO-DIMENSIONAL FLOW 


The basic difficulty in the two-dimensional case is due to the large discretization 
error arising from an inconsistent treatment of shocks and other discontiritities. 
This error will remain large as the mesh is refined except in a few special cases 
which are discussed below in a sinplified manner. For this case the system of 
conservation laws is 


9 F ^ „ 

9t ax ay 


( 24 ) 


where F= { p, m, n, e} , H = { m, — +p, 

2 

G= {n, is the y-momentum, and 


(e»p) - } , 


P= (y-l) (e - j (m2 +n2) /p) . 
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If shocks appear, they satisfy the steady state jimp condition 

^ [H]^ + tG]y = 0 (25) 

where dy /dx is the local shock inclination and are the jmps 

^ X 

in the x , y - directions . 

The MC scheme for eq. (24) is: 







a h (G.”. -g”, ) 

yy'-j,k j.k-oJ 


(26 a) 






a h 

X X 


< 1 . 



n+1 n+1 


(26 b) 


where h = At /Ax, h = At /Ay, and o , a =±1 determines the differencing 
X y X y 

sequence in the two directions. 


The error modes are excited by the discretization errors T. , (eqs. (6) and (7)) 

J > 

which act as local source terms. In the scheme considered the magnitude of these 

errors are 0( Ax^, Ay^) and can thus be made arbitrarily small in smooth regions 

of the field be refining the mesh. This is, however, not the case for grid points 

close to discontinuities, and in the following we shall only consider the part 

of T. , which is not necessarily 0( Ax^, Ay^). For the MC scheme this part is 
J » ^ 


T. 


j,k " ■ " %-l,k ^ V ^%,k+l "%j,k-l^ ^ 




(27) 


where h = Ax /Ay, subscriot 0 means the exact steady solution as before, 

and aHq, a Gq is the change in H, G after application of the predictor step 

on the exact solution. If the change in the Jacobians A. . and B. , caused by 

J > 1 

the predictor step is neglected, we have: 
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"^j.k " ^)j,k ‘^j.k 

Boj.k^j.k 

®j.k * ”'='x^c%,k ““Oj-a^.k^ ■ °yV^^j.k‘^0j,k-Oy) 

For sin|)licity the case of two constant parallel flow regions separated by a 
straight shock will be considered. The local difference between this case and 
a smoothly varying field is only due to terms vdiich vanish as the mesh is refined. 
The term in square brackets in eq. (27) for T- , can only vanish in two special 
cases. If the mesh is aligned with the shocks the junps in H and G across the 
shock are zero. For this case the ranaining terms also vanish, and the previous 
oie-dimensiorial results are applicable. 

Secondly the first term of eq. (27) is also ccarsistent with the shock junp rela- 
tion (25) if the local inclination of the shock <fy^/dx is equal to Vh^ , in 
vdrich case the shock is parallel to the mesh diagonals. To make the remaining 
terms in eq. (27) vanish the differencing sequence cannot be arbitrarily chosen, 

i.e. a , a . From eq. (28) it follows that differences should be taken across 

y 

the shock in both directiois in either the predictor or the corrector step, and 
not in one direction in each step. 

If neither of these two ccaiditions are satisfied, the discretization error is 
always 0 (1) regardless of mesh refinements. The largest errors are produced 
when the scheme captures the jump in only one direction. This will occur at 
regular inteivals along a straight, oblique shock in a rectangular mesh. 

A good exanple of these effects is seen in the ccnputations of Kutler et. al. [10] 
for the reflection of an incoming shock on a wedge in supersonic flow. The errors 
are large along all shocks except for the reflection from the wedge surface, where 
the above conditions seem to be at least approximately satisfied. The results 
given here are easily extended to self-similar problems such as [10] . Another 
exanple for a self-similar case is also given in [11] for the reflection of a 
strong shock from a wedge. The mesh and differencing sequence was chosen to mini- 
mize the errors along the discontinuities and the result is shown in Fig. 6 . 
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For the MC scheme with permuted differencing [ 5 ] , the time-split scheme 1 8 ] and 
the iqjwind WB schane the square bracketed term of eq. (27) is the same as for the 
regular MC scheme. The remaining terms are, however, always 0(1) for these 
schanes. It is therefore possible that the advantage of these schemes resulting 
from better stability and faster convergence towards the steady state, may be 
largely offset by the better shock-capturing properties of the regular MC scheme 
vdien properly applied. 


5. CONCLUDING REMARKS 


It has been explained in the previous sections how the discretization errors near 
steady shocks can excite the spurious e3qx>nential error modes of a shock- capturing 
scheme, and also vdiy these errors can exist even in the special cases when the 
discretization errors are small. We feel that proper understanding of the be- 
haviour of the spurious modes is inportant to improve the quality of the steady 
state results obtained from such schemes. 

In the arguments above it has been assumed that the shock is correctly located in 
the mesh. A wrong shock location can be interpreted as a large discretization 
error locally and will hence lead to large errors in the vicinity of the shock. 

To minimize this possibility a fully caiservative scheme should be used. It is 
worth noting, however, that a conservative scheme does not guarantee correct shock 
locations. According to the theory above a wrong shock location in the steady 
state is also possible with the use of a shock-fitting technique, provided that 
the error introduced can decay sufficiently fast towards the boundaries. The reason 
why the present two-dimensional shock- fitting techniques [12] work so well is 
probably because the special numerical schemes used along the shock boundary tend 
to suppress the oscillating exponential error modes [13] . The applicability of 
such schemes is unfortunately severely limited for three-dimensional cases. 

The shock-capturing approach with conservative schemes can be applied without 
difficulty in problems of arbitrary dimension. The results given here are also 
easily extended to multi-dimensional cases. 
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Table I. Comparison of relative errors 


j 

Mass flux error 
•10^ 

Theory Conrouted 

Momentum flux error 
•10^ 

Theory Computed 

Energy flux error 
•10^ 

Theory Computed 

23 

-4.485 

-4.440 

-1.111 

-1 .079 

-2.701 

-2.647 

24 

6.003 

5.898 

-1 .003 

1 .045 

3.814 

3.740 

25 

-5.912 

-5.927 

1.132 

1 .098 

-4. 463 

-4. 449 

26 

4.505 

4.505 

-0.862 

-0.886 

3.401 

3.417 

27 

-3.433 

-3.433 

0.657 

0.639 

-2.592 

-2.577 

28 

2.616 

2.598 

-0.501 

-0.514 

1 .975 

1 .976 
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Fig. 6 Reflection of shock wave from a wedge. Density contours. 
Full lines are drawn from a Schlieren photograph of the 
corresponding shock tube experiment. 


